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A three-dimensional granular system fluidized by vertical container vibrations was studied using 
pulsed field gradient (PFG) NMR coupled with one-dimensional magnetic resonance imaging (MRI) . 
The system consisted of mustard seeds vibrated vertically at 50 Hz, and the number of layers Nt < 4 
was sufficiently low to achieve a nearly time-independent granular fluid. Using NMR, the vertical 
profiles of density and granular temperature were directly measured, along with the distributions of 
vertical and horizontal grain velocities. The velocity distributions showed modest deviations from 
Maxwell-Boltzmann statistics, except for the vertical velocity distribution near the sample bottom 
which was highly skewed and non-Gaussian. Data taken for three values of Ni and two dimensionless 
accelerations V = 15, 18 were fit to a hydrodynamic theory, which successfully models the density 
and temperature profiles including a temperature inversion near the free upper surface. 

PACS numbers: 45.70.Mg, 76.60.Pc, 81.05.Rm 



I. INTRODUCTION 

One of the basic granular flow phenomena is the cre- 
ation of a fluidized state by vibration of the container 
that holds the granular medium. If the container is 
vibrated vertically with a sinusoidal waveform z(t) — 
ZQCOs(ujt), the dimensionless acceleration amplitude rel- 
ative to gravity is T = uj 2 zo/g and T > 1 is required to 
set the granular system into motion 1]. A second key 
dimensionless parameter is X — Ne(l — e), where e is 
the velocity restitution coefficient for grain-grain colli- 
sions and Ni is the number of grain layers (the height 
of the granular bed at rest divided by the grain diame- 
ter). Experiments and simulations in one and two dimen- 
sions suggest that X < Xf with Xf ~ 1 is required to 
reach a uniformly fluidized state, in which the stochas- 
tic motions of the grains exceed motions that are har- 
monically or sub-harmonically related to the container 
vibration [|- [3-0 For larger X > Xf inelastic collisions 
rapidly quench stochastic grain motion and the granular 
bed moves as a coherent mass over much of the vibration 
cycle, leading to phenomena such as subharmonic surface 
patterns and localized "oscillon" structures H, 0, 13 • 

In this paper we report an experimental study of the 
uniformly fluidized state in three dimensions, using NMR 
techniques to probe the grain density profile and velocity 
distribution. Many authors have applied the methods of 
statistical mechanics to the fluidized granular state |9LIT(i 
ITU Il2l IT3L IT3. ITU lill . Typically, one defines a "granular 
temperature" proportional to the random kinetic energy 
per degree of freedom of the grains, and hopes to apply 
hydrodynamic concepts such as thermal conductivity and 
viscosity to the granular fluid. Using NMR we are able 
to directly measure the granular temperature profile in 
a dense, three-dimensional sample for comparison with 
hydrodynamic theory. 

Such comparisons are significant as the applicability of 



hydrodynamics to granular systems is questionable. In 
a vibrofluidized state with energy injection from below, 
as studied here, the granular temperature varies sharply 
with height over distances comparable to the mean free 
path. Hence the separation between microscopic and 
macroscopic (hydrodynamic) length scales is marginal, 
and hydrodynamics may or may not provide an accurate 
description of the system 0, 0, . 

Previous experiments have created (quasi) two- 
dimensional systems by confining grains between verti- 
cal plates or on a horizontal or tilted surface, enab ling 
data collection by high-speed video [H |H HI HI 
I2EL I26L |27| . Three-dimensional granular systems have 
been studied experimentally by external probes such as 
capacitance [28|], mutual inductance [2ij and diffusing- 
wave spectroscopy [3(|, as well as by non- invasive inter- 
nal probes including NMR [U 0Jp and 
positron-emission particle tracking (PEPT) 37]. The 
PEPT method has been used to measure the granular 
temperature profile of a vibrofluidized bed at volume fill- 
ing fractions up to 0.15 37]. In the present study using 
NMR, we have measured granular temperature profiles 
at much higher densities, up to 0.45 volume filling frac- 
tion. Density is a key parameter for granular flows, as 
studies suggest the existence of long-lived fluctuations 
and glassy, out-of-equilibrium behavior as the density in- 
creases towards the random-close-packed volume filling 
fraction rs 0.6, si g nalin g a breakdown of conventional 
hydrodynamics |38|, l39l |4(J, |4l|. Even in the absence 
of glassiness there are precollision velocity correlations 
which increase with density [Til . E^ : these correlations 
are included in some hydrodynamic theories |13| but ig- 
nored in others ("molecular chaos" assumption) |l5j . 

Other workers have used NMR to study the dense 
granular flows that occur in a rotating drum [M], [3(|, 
for quasi-static tapping [32), in a period-doubled arch- 
ing flow and in a Couette shear apparatus j3f|- In 
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contrast to these earlier studies we probe a system suf- 
ficiently simple to be described by first-principles gran- 
ular hydrodynamics. In addition, we use rapid gradient 
pulse sequences to probe deep within the ballistic regime, 
which has not been possible for these very dense systems. 

In the experiments described here, small samples con- 
sisting of 45 - 80 mustard seeds of mean diameter 1.84 
mm were confined to a cylindrical sample tube with 9 mm 
inside diameter, and vibrated vertically at oj/2-k = 50 Hz, 
r = 15-18. The sample container was sufficiently tall 
to prevent grain collisions with the container top. The 
total number of grains is small in our experiments (com- 
pared, for example, with the number of particles in recent 
molecular-dynamics simulations |43j). However these 
physical experiments include effects present in practical 
applications but not usually addressed by simulations: 
irregular and variable grain size, grain rotation, realistic 
inelastic collisions, and the effect of interstitial gas. 

The granular systems studied here were small both hor- 
izontally (5 grain diameters) and vertically (Ng < 4), and 
thus it may be surprising that they can be successfully 
described by a hydrodynamic theory. The small sample 
diameter was simply a consequence of the inner diameter 
of the NMR apparatus that was available, and could be 
increased in future studies. Conversely, the sample depth 
is inherently limited to a small number of layers by the 
fluidization condition Ng(l — e) < 1. Therefore physi- 
cal vibrofluidized systems (which typically have e < 0.9) 
are limited to Ng < 10 and the applicability of granular 
hydrodynamics is a nontrivial question to be tested by 
experiment. 

A short report of some of this work was published in 
Ref. l44l using a single set of Ng, T values. For the present 
report, both Ng and T were varied. In addition improved 
methods were developed for calibrating the NMR gradi- 
ent coils, resulting in more accurate granular temperature 
profiles. In particular, better agreement is found with hy- 
drodynamic predictions for the "temperature inversion" 
that occurs at the free upper surface. 



II. EXPERIMENTAL METHODS 

A. Experimental setup for NMR of vibrated media 

The experiments were carried out in a vertical-bore 
superconducting magnet with a lab-built 40 MHz NMR 
probe equipped with gradient coils for pulsed field gradi- 
ent (PFG) and magnetic resonance imaging (MRI) stud- 
ies. The low static field value (0.94 T) was chosen to min- 
imize magnetic field gradients in the sample due to the 
susceptibility difference between grains and voids. This 
permitted the use of a simple unipolar PFG sequence 
(Fig. nj for the fastest possible time resolution. The 
9 mm inside diameter, 6.5 cm high sample chamber had 
a glass side wall to minimize energy transfers between 
grains and the side wall. The teflon bottom wall was 
pierced by an array of 30 0.46 mm diameter holes, making 



it possible to evacuate the sample container. An exter- 
nal vacuum gauge registered less than 200 mtorr during 
experiments with evacuated samples. 

A loudspeaker connected to the chamber by a fiberglass 
support tube and driven by a function generator and au- 
dio amplifier was used to vibrate the sample container 
vertically with a sinusoidal waveform. The amplitude 
and phase of the container motion were precisely mea- 
sured by fitting the digitized output of a micromachincd 
accelerometer (Analog Devices ADXL50) mounted to the 
support tube just below the NMR magnet. The gradient 
coils used for position (MRI) and displacement (PFG) 
measurements were calibrated by combining data from 
measurements on phantoms of known dimensions, along 
with diffusion measurements on a water sample. In this 
way we obtained calibration maps for each gradient coil 
of both the gradient in the designed direction, and the 
total gradient magnitude including components orthog- 
onal to the designed direction. The accelerometer and 
vertical gradient calibrations were cross-checked by mea- 
suring the displacement waveform for a granular sample 
consolidated with epoxy. The dimensionless acceleration 
values r reported here are accurate to within ±0.1. 



B. Granular samples 

Previous NMR studies of granular media used plant 
seeds or liquid-filled capsules to take advantage of the 
lon g-liv ed NMR signals produced by liquids 0, |U IH, 
133. l35l |36l. For the present study, we chose mustard 
seeds [45] which showed no tendency to static charging. 
The seeds were ellipsoidal, with three unequal major di- 
ameters and eccentricity < 20%. The average of all three 
diameters measured for a sample of 10 seeds was d = 1.84 
mm with a standard deviation of 5%, and the average 
mass was m — 4.98 mg/seed. These average values of 
d and m are used below for comparison with a hydrody- 
namic theory, which is formulated for samples of identical 
grains. 

To roughly measure the restitution coefficient e, a video 
camera was used to record the bounce height when the 
seeds were dropped onto a stone surface. As some linear 
kinetic energy is converted to rotational energy upon im- 
pact, this gives a lower bound on e. Conversely the small- 
est drop velocity that could readily be used (44 cm/s) was 
larger than the typical RMS velocity in our vibrofluidized 
experiments (< 30 cm/s). Since e is often a decreasing 
function of impact velocity [4(|[47|], the restitution coeffi- 
cient in the vibrofluidized experiments should be greater 
than is measured in the dropping trials. The e value de- 
duced from the highest rebound over 15 dropping trials 
on each of five seeds ranged from 0.83 to 0.92, with an 
average of 0.86. 
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can be repeated as many times as necessary to build up 
good signal/noise. Allowing full longitudinal relaxation 
between acquisitions (the simplest but not the fastest way 
to take data), the acquisition time for a time-resolved 
propagator p(Au, v, t) was 21 hours. 



D. Measurement of RMS displacement 



FIG. 1: Pulse sequence used for PFG/MRI data acquisi- 
tion. This is a stimulated-echo PFG sequence followed by one- 
dimensional frequency-encoded MRI. The three solid vertical 
bars are ^ RF pulses, the boxes with horizontal hatching are 
gradient pulses in the displacement-resolution (PFG) direc- 
tion, and the open boxes are gradient pulses in the position- 
resolution (MRI) direction. The gradient pulses are labeled 
with the gradient strengths in T/m, while numbers on the 
time axis indicate durations of intervals in microseconds. The 
start of the sequence was triggered at a fixed phase of the sam- 
ple vibration, and the delay time Dl was varied over the range 
1-20 ms to select the time in the vibration cycle at which NMR 
data was acquired. The displacement measurement time At 
is equal to the time between PFG gradient pulses plus one 
third the duration of the gradient pulses; for the sequence 
shown At = D2 + 380 /js. An eight-sequence phase cycle was 
used to suppress undesired coherences and to enable acqui- 
sition of both real and imaginary parts of the complex PFG 
propagator. 



C. Data acquisition 

Proton NMR signals from the mustard-seed sample 
were acquired using the PFG/MRI sequence shown in 
Fig. ^ After Fourier-transformation with respect to 
time, this sequence yields the position-dependent PFG 
signal (one-dimensional image) 



h(q u ,v) = / p(Au,v)e l ^ u dq u 



(1) 



where v — x, y or z is the direction of the MRI gradient 
and u = x, y or z is the direction of the PFG gradi- 
ent. In this equation q u is proportional to the area of the 
PFG gradient pulse, which is varied through an array of 
regularly spaced values (Fig. . Here the "propagator" 
p(Au, v) is the joint probability that a grain (more pre- 
cisely an oil molecule in a grain) is located at position 
v and moves a distance Au during the time interval At 
set by the pulse sequence • From Eq. ^ the propaga- 
tor p(Au, v) can be obtained by Fourier transforming the 
NMR signal h(q u ,v) with respect to q u . With the gra- 
dient values used in these experiments the propagator 
was obtained with resolutions of 500 /mi in the position 
variable v and 42 [im in the displacement variable u. 

This type of NMR experiment does not determine the 
positions or displacements of individual grains. Rather, 
an ensemble average, the propagator, is obtained. Pro- 
vided that the propagator has a well defined time average 
at each phase of the vibration cycle, the NMR acquisition 



For some purposes, the full propagator p(Au, v) is of 
interest, while for other purposes (e.g. measuring the 
granular temperature) only the position-dependent RMS 
displacement is required, 
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(2) 



Here n(v) — J p(Au, v)d(Au) is the position dependent 
number density. Both n(v) and Au rms (v) can be ob- 
tained directly from the NMR data (before it is Fourier 
transformed to yield the propagator) using 
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For optimum signal/noise we evaluated the RHS of these 
expressions by fitting the central (small \q u \) portion of 
h(q u ,v) to a function including constant, quadratic, and 
higher terms. 



III. EXPERIMENTAL RESULTS: DENSITY 
AND DISPLACEMENT DISTRIBUTIONS 

In the present experiments, all measured quantities de- 
pended strongly on height z, due to gravity and the local- 
ization of energy injection at the sample bottom. There- 
fore, most data were taken with the MRI axis v — z. As 
both vertical and horizontal grain motion were of inter- 
est, for the PFG axis both u = z and u — x were used. 
Thus, the primary experimental data shown in this paper 
are the height-dependent number density n(z), the ver- 
tical and horizontal displacement propagators p(Az, z) 
and p(Ax,z), and the vertical and horizontal RMS dis- 
placements Az Tms (z) and Ai rms (z). In addition, some 
data were taken using a two-dimensional MRI sequence 
to verify that the grain density was approximately uni- 
form in the horizontal plane and to calibrate the number 
of layers Nt (Appendix . 

The granular state with Nt = 3.0, T = 15, and sample 
cell evacuated was selected as a reference state measured 
in greatest detail, and other states were selected by vary- 
ing a single parameter from the reference value. The 
following variations were carried out: (a) increase of am- 
bient air pressure up to 1 atm, (b) variation of Nt from 
2.2 to 4.0, (c) variation of T up to 18. The upper limit on 
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FIG. 2: Effect of cell evacuation. Un-normalized number 
density (top) and mean-square vertical displacement in time 
interval At = 1.38 ms (bottom) as functions of height z for 
Ni — 3.0 layers of mustard seeds vibrated at lo/2h — 50 Hz 
and acceleration amplitude relative to gravity T = 15. The 
data labeled Air were taken in normal atmospheric condi- 
tions, while for the data labeled Vacuum the cell was evacu- 
ated so that an external pressure gauge registered less than 
200 mtorr. The small differences between the Air and Vacuum 
data are comparable to the typical experimental repeatability 
and therefore are not significant. 



r was set by the apparatus, while the lower limit on T and 
the upper limit on N# were set by the requirement that 
the sample fully and consistently fluidize, when observed 
visually outside of the NMR apparatus. For smaller T 
or larger Ne we observe a complex, hysteretic transition 
to a partially fluidized, partially jammed state which is 
beyond the scope of the present study. 



A. Effects of Interstitial Air 

Most of our data were taken with the sample container 
evacuated. Some of the measurements were repeated in 
normal atmospheric conditions, to determine how large 
the interstitial gas effects could be. Experimentally, the 
effect of air is unobservably small (Fig. |2J). This agrees 
with the following rough estimate of viscous energy losses 
due to air drag. The mean free path varies widely with 
height in our system but a typical value is A ~ 0.3d 
(Appendix lB|l . Similarly, a typical grain velocity is on 




n (a.u.) 



FIG. 3: Un-normalized number density of grains n as a func- 
tion of height z and time in the vibration cycle t for the ref- 
erence vibrofhhdized state Ne = 3.0, F = 15. The thick black 
curve shows the container vibration waveform, obtained by 
fitting the digitized accelerometer output. For this figure the 
curve is drawn to show the motion of the container bottom, 
which impacts the sample bottom at t ~ 12 ms. Apart from 
the disturbance caused by this impact, the density profile is 
largely independent of t. 



the order of 20 cm/s. Using these values and the drag 
force on a sphere in an unbounded medium, we compute 
that the fraction of kinetic energy lost to viscous drag 
over one mean free path is 9 x 10~ 4 for p — 1 atm and 
3 x 10~ 4 for p = 200 mtorr. These energy losses are 
negligible compared to the fraction of kinetic energy lost 
in an intergrain collision, 1 — e 2 w 0.24. 



B. Time-resolved density and displacement profiles 

In Sec. II VI below the NMR data are fitted to a hydro- 
dynamic theory which models the granular system as a 
fluid with time-independent (but spatially varying) grain 
velocity distribution. In this picture, the vibrating con- 
tainer bottom serves to inject energy (granular "heat") 
into the system. For such a picture to make sense, the 
overall state of the granular system must be nearly inde- 
pendent of time within the vibration cycle. 

We have taken some NMR data as a function of t 
within the 20 ms container vibration period, Figs. and 
0] The density profile (Fig. |3Jl is indeed largely indepen- 
dent of t over most of the sample, except for a noticeable 
disturbance at small z for t near the top of the container's 
motion, when the container floor impacts the sample. In 
the bottom layer of the sample, this disturbance damps 



FIG. 4: Profiles of the vertical (left) and horizontal (right) RMS grain displacement Az rms , Az rms in a time interval At = 
1.38 ms, as functions of the height z and the time in the vibration cycle t. The experimental conditions are the same as for 

Fig. on 



out (thermalizes) over a time span of less than 10 ms, and 
above the bottom layer it is hardly evident at all. Sim- 
ilarly, in the RMS displacement profile (Fig. 0} a large 
transient is visible in the vertical displacement due to 
grains impacted by the container bottom, while the hori- 
zontal displacement transient is much smaller relative to 
the time-average value. 

Over most of the cycle the sample "floats" well above 
the vibrating container floor (Fig. EJ) , illustrating two im- 
portant points about the vibrofluidization process: (a) It 
is possible to achieve a nearly time-independent fluid- 
like granular state even though the driving vibration 
amplitude is much greater than typical stochastic mo- 
tion within the sample. The required conditions are 
(i) uj* — ^(d/g) 1 / 2 > 1, which ensures that the sam- 
ple itself only falls a grain-scale distance in one vibration 
period no matter how large the container oscillation am- 
plitude is, and (ii) X = Nt(l — e) < 1, allowing the sam- 
ple to fully fluidize. For the experiments described here 
uj* = 4.3 and 0.29 < X < 0.52. (b) For realistic multi- 
layer granular systems, the primary contact between the 
grains and the vibrating cell bottom occurs near the top 
of the vibration stoke, when the cell bottom velocity is 
much less than its RMS value. Therefore, the energy 
transfer into the granular medium cannot be estimated, 
even roughly, from the RMS driver velocity. This is be- 
cause the energy transfer depends upon the mean cell 
bottom velocity during grain impacts |49l |. 

Apart from Figs. and 0] the data shown in this paper 
have all been averaged over the vibration cycle for im- 
proved signal/noise. It should be born in mind that the 
state of the system is not, in fact time-independent very 
near the sample bottom and that the complex process of 



energy transfer from the vibrating base is not captured 
by these time averages. 



C. Short time displacements: Crossover to the 
ballistic regime 

Figure El shows the RMS vertical and horizontal grain 
motion in the reference state Ni = 3.0, V = 15 as func- 
tions of the displacement measurement time At, for slices 
at three different heights z within the sample. At each 
height there is a crossover at short At to the power law 
Az rms , Aa; rms tx At, which indicates ballistic behavior 
(displacements proportional to time). For longer At the 
data fall below the ballistic power law, possibly crossing 
over to a diffusive power law Az rms , Aa; rms cx (At) 1 / 2 . 
The crossover time, which is an estimate of a typical 
grain collision time, ranges from 5 ms in the dense, lower 
portion of the sample to 10 ms in the less dense, up- 
per portion of the sample. In Appendix [5] it is shown 
that the crossover time agrees well with the grain col- 
lision time computed from the mean free path and the 
measured granular temperature. 

A key feature of our experiments is the ability to take 
data for displacement intervals At much less than the 
grain collision time, i.e. deep within the ballistic regime. 
The RMS displacement data shown in the remainder of 
this paper was taken with At = 1.38 ms, well within the 
ballistic regime for all values of Ng,F used. 
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At (ms) 

FIG. 5: Root-mean-square vertical (circles) and horizontal (squares) grain displacements Az rms , Ax rms versus the displacement 
measurement time At for a sample with Ne = 3.0, V — 15. The data are shown for three different horizontal slices (z ranges), 
which may be correlated with the density profile in Fig. |3 In each log- log plot, the solid line has a slope of unity (RMS 
displacement proportional to time) while the dashed line has a slope of |. 



D. Horizontal and vertical velocity distributions: 
Deviations from Gaussian 

Figures El and [7| show the distributions of vertical and 
horizontal displacements, respectively, for the reference 
vibrofluidized state Ng = 3.0, T = 15 and displace- 
ment observation time At = 1.38 ms within the ballistic 
regime. These propagators should be good approxima- 
tions to the distributions of vertical and horizontal ve- 
locities, p(v ZjX ,z) with v z , x = Az, xj At. 

In some previous experiments 12H. I25I |26| and sim- 
ulations [50| of granular fluids (mostly in reduced dimen- 
sions) significant deviations have been found between 
the measured velocity distribution, and the Gaussian 
(Maxwell-Boltzmann) distribution that describes an elas- 
tic fluid in equilibrium, p(vi) cx exp(— mvf/2T) (setting 
Boltzmann's constant to unity). For example, Rouyer 
and Menon [2{| found that the distribution of horizon- 
tal velocities in a two-dimensional vibrofluidized media 
robustly obeyed 

p(v x ) (X ex.p[-\v x /a a \ a ] (4) 

with a = 1.55 ± 0.1. Using kinetic theory it has been 
shown that the tails of the velocity distribution should 
follow Eq.^Jwith a = 1 for a freely-cooling granular gas, 
and a = 3/2 for a granular gas with stochastic forcing ap- 
plied to every grain |51|. However, for the freely-cooling 
case the central part of the distribution is nearly Gaus- 
sian (a = 2) and the crossover to a = 1 occurs very far 
out in the tails [52|. 

In the present experiments and those of Ref. |2{| a 
steady state was achieved by supplying energy from the 



boundaries of the system, which corresponds neither to 
the freely cooling nor to the uniformly heated case. A 
few recent theoretical and simulation studies have treated 
the two-dimensional granular fluid with a strong thermal 
gradient, as occurs in vibrofluidization experiments, and 
have found features such as anisotropic granular tem- 
perature, asymmetric velocity distributions, and non- 
universal values of the velocity-stretching exponent a 

Usui mm. 

In previous work the measured velocity distribution 
typically has been scaled by the RMS velocity a (note 
a 7^ a a)- This is possible for simulations and for video 
experiments in which every particle is tracked, allowing 
a to be computed unambiguously. In NMR experiments 
there is a noise floor, visible in Figs. 16171 which could con- 
ceal tails of the distribution contributing to a. One pos- 
sibility would be to use direct measurement of the RMS 
velocity from q-space data, described in Sect. Ill Dl and 
used below to measure the granular temperature. Here 
we have taken the simpler approach of directly fitting the 
experimental velocity distributions to Eq. 0] allowing the 
width parameter a a to vary. The results of these fits are 
shown in Fig. [5] 

The fitted stretching exponents for vertical and hori- 
zontal velocity distributions a z , a x generally lie between 
1.5 and 2.0, with values closer to 2.0 (Gaussian) near 
the top of the sample. Altough a z and a x are not sig- 
nificantly different, the large error bars for a z near the 
sample bottom reflect the fact that the experimental ve- 
locity distribution is asymmetric (Fig. [HJ and hence does 
not fit Eq. 0] well. The present results for a x in a three- 
dimensional vibrofluidized system appear marginally dif- 
ferent from the corresponding two-dimensional results of 
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-0.05 0.05 

FIG. 6: Propagator, or distribution p(Az, z) of vertical displacements Az during time interval At — 1.38 ms, as a function of 
height z for the reference vibrofluidized state Ne = 3.0, T = 15. As At is well within the ballistic regime, p(Az, z) approximates 
the distribution of vertical velocities p(v z , z) with v z = Az/ At. The same data are shown from two different viewpoints in the 
left and right plots. Near the bottom of the sample (small z, most visible in the left plot) the velocity distribution is highly 
skewed, with a wide tail on the v z > side due to particles traveling upward with large velocity after impacting the vibrating 
container bottom. Well away from the bottom (right plot) the velocity distribution is nearly symmetric. 




-0.05 0.05 

FIG. 7: Two views of the propagator p(Ax, z) for horizontal displacements Ax in the same conditions as in Fig. |S| The 
distribution of horizontal velocities is symmetric at all heights z, consistent with the symmetry of the experimental arrangement. 



Ref. |25[ Some of this difference may simply be due to the 
different way in which we fit the data, allowing a a to vary 
rather than fixing it to a precisely measured experimen- 
tal value. Our results arc similar to the two-dimensional 
simulation results of Ref. |50| . which found deviations 
from Gaussian are most pronounced near the bottom of 
the system. 



IV. FIT TO A HYDRODYNAMIC THEORY 

A. Granular hydrodynamics for a stationary 
vibrofluidized state 

Granular hydrodynamic equations have been derived 
by many groups, using increasingly sophisticated approx- 
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FIG. 8: Fitted stretching exponents a z<x for the verti- 
cal and horizontal velocity distributions versus height z, for 
Nt = 3.0, r = 15. The Gaussian and exp(— v 3 ^ 2 ) cases are 
shown by the dashed lines. The large uncertainties for a z at 
small z indicate a poor fit of the stretched-exponential func- 
tion, which is symmetric, to the highly asymmetric velocity 
distribution (Fig.|BJ. 

imation schemes [1 H 111 III III III El In this section 
we show how any of these theories may be numerically in- 
tegrated for comparison with experiments. Section llV Bl 
shows this comparison between the theory of Ref . Il5l and 
our experimental results. 

In Ref. |4jj | granular hydrodynamics was successfully 
used to describe a nonstationary simulation with large 
variations of the density profile over the course of a vi- 
bration cycle. Here, we specialize to the simpler situation 
in which the average velocity field v(r, t) vanishes and the 
number density field n(r, t) and temperature field T(r, t) 
are independent of t and depend upon r only via the 
height coordinate z. We consider a system of N identical 
spherical grains of diameter d and mass m. 

With these assumptions the pressure p{z) and number 
density n(z) satisfy a hydrostatic equilibrium equation 

dp/dz = —mgn (5) 

and the upward heat flux q z (z) obeys an energy flux bal- 
ance equation 

dq z /dz = —P c (6) 

where P c is rate of kinetic energy loss per unit volume, 
due to inelastic collisions. The heat flux is 

q z = — n(dT j dz) — /i(dn/dz) (7) 

where n is the conventional thermal conductivity, and /i 
is a transport coefficient that is nonzero only for inelas- 
tic systems ^2- Due to the second term in Eq. \7\ it 
is possible to have a flow of heat from colder regions to 



warmer ones. Although a number of authors have cal- 
culated a [Til Ibl llH and discussed the consequences of 
/i ^ |5rl l57l l5Sj|. to our knowledge no heuristic ex- 
planation has been proposed for the existence of a heat 
current driven by density gradients. This density-driven 
heat current has a marked effect in the present experi- 
ments. 

Equations 15171 are simply the results of definitions and 
of momentum and energy conservation. A hydrodynamic 
theory closes this set of equations by providing expres- 
sions for p(n,T) (the equation of state), P c , k and \i. 
It is convenient to nondimcnsionalize variables using the 
grain diameter d, grain mass m, and gravity g, 

z* = {z-z )/d, (8) 

ri* = nd 3 , 

T (*) = T {i)/mgd, 

q* z = q z {d/gfl 2 /m. 

Here zq is the height of the bottom of the sample, corre- 
sponding to z* = 0. 

For a hard-sphere system with constant restitution co- 
efficient the equation of state takes the form 

p/nT = f(n*). (9) 

Using Eas. 15171 and l9l the following expression is derived 
for the density gradient: 

dn = cIz/k - mg/f 

dz ( T /n)(l + f&)-iM/K 

Equations |5l |S] and ^] are coupled first order ordi- 
nary differential equations for {p, q z ,n} that can be nu- 
merically integrated from the sample bottom at zq to 
z = +oo. The boundary conditions are as follows. For 
a system in gravity with no upper wall, the pressure and 
heat current satisfy p(+oo) = 0, q z (+oo) = 0. Using 

Eq. [SJthis gives p(zo) — rngriA, where tia = n(z)dz 
is the total number of grains in the sample per unit hor- 
izontal cross section. The heat current into the system 
at the bottom q z (zo) is an input parameter determined 
by the vibration amplitude and frequency and the sys- 
tem response. For given values of q z (zo) and p(zq), the 
density at the bottom n(zo) must be adjusted to satisfy 
the upper boundary condition q z (+oo) = 0. 

Using our convention, the dimcnsionless (starred) 
forms of Eqs. 151 1UI and the boundary conditions are ob- 
tained by setting m = g = d = 1. We have carried out 
the numerical integrations and present comparisons with 
the experiments in this dimensionless form. 

B. Fit of the experimental results to the 
Garzo-Dufty hydrodynamic theory 

In this section, we show fits of our experimental data to 
a formulation of granular hydrodynamics by Garzo and 
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FIG. 9: Experimental density and granular temperature pro- 
files (symbols) and fit to the hydrodynamic theory of Ref. |l5| 
(curves) for Nt = 3.0 layers of grains vibrated at dimension- 
less acceleration amplitude F = 15. Dimensionless variables 
defined in Eq.|H]are used. In the top graph the symbols labeled 
n* show the density profile n*(z*) derived from the data set 
used to measure vertical displacements while the symbols la- 
beled n% show the same quantity derived from the horizontal- 
displacement data set; these two profiles should be identical 
and differences reflect experimental accuracy. In the bottom 
graph T*, T* are the measured granular temperature profiles 
for vertical and horizontal displacements, respectively, mea- 
sured using observation interval At — 1.38 ms. It is expected 
that T* > T* , at least near the bottom of the sample, as 
energy is input only to the vertical degrees of freedom by the 
contain vibration. Fhe theory should be compared with the 
symbols labeled T* , which are computed from the measured 
temperatures using T* = gTj + §T£. 

Dufty . This calculation is based on the "revised En- 
skog theory" , which is accurate for elastic fluids over a 
wide range of densities and is extended in Ref. to ar- 
bitrary values of the restitution coefficient e. Some other 
theories of granular hydrodynamics are derived as expan- 
sions in powers of 1 — e 2 [lj( , and thus may only be accu- 
rate for £ 1. On the other hand, the Garzo-Dufty the- 
ory makes the molecular chaos assumption that the pre- 
collision grain velocities are uncorrelated. Pre-collision 
velocity correlations are known to exist in inelastic hard- 
sphere fluids 0, and can be treated theoretically 
using "ring kinetic theory" [lflj . 

It is doubtful whether fits to the present data could 



FIG. 10: Density and temperature profiles for Ne = 2.2, with 
other conditions the same as in Fig. [5] Under these conditions, 
a significant discrepancy can be seen between the experimen- 
tal and theoretical granular temperatures in the lower portion 
of the sample, z* < 3. We observe this discrepancy in condi- 
tions such that the density in the lower portion of the sample 
is well below the peak density (upper graph) . 



distinguish between these various hydrodynamic theories. 
We chose to compare our data to Ref. because this cal- 
culation is nominally valid for wide ranges of density and 
inelasticity, and because the theoretical results for f,P c , 
k and /i are presented in an algebraic form convenient 
for numerical integration. In this theory the granular 
temperature is assumed to be isotropic, T x = T y = T z , 
where Ti = m(vf), but it is not assumed that the velocity 
distribution is Gaussian. Experimentally we find modest 
differences (of order 20% or less) between the horizon- 
tal and vertical granular temperatures, T x ^ v ^ T z , and 
more generalized formulations of hydrodynamics may be 
able to treat this anisotropy |53l |5J, |53 • For comparison 
with the isotropic-temperature theory the measured ver- 
tical temperature T z and horizontal temperature T x are 
averaged using T — |T Z + |T X . 

Figures El through ^] show a joint fit to the Garzo- 
Dufty theory of our data for three different values of the 
dimensionless bed height N( and two different values of 
r. For this fit the height of the sample bottom zq was 
allowed to be different for each data set, since the sample 
"floats" above the vibration center by an amount that de- 
pends upon the dynamic interplay between energy input 



10 





0.8 
0.6 
* s 0.4 
0.2 

0.0 

2 



* 1 





-2 2 4 6 8 10 

z* 

FIG. 11: Density and temperature profiles for Ne = 4.0, 
with other conditions the same as in Fig. [3] The peak density 
n* = 0.85 observed in these conditions corresponds to volume- 
filling fraction cj> = (n/6)n* — 0.45, about 75% of the random 
close packed value. 

and gravity (Fig. |3J). Similarly, the energy input ^(0) 
was a separate fitting parameter for each data set. Con- 
versely a single value of the restitution coefficient e was 
used to fit all four data sets. Thus the only fitting pa- 
rameters are zq, q*{fS) and e. All other quantities entering 
the fits shown in Figs. 1911^1 are measured experimentally, 
including all of the axis variables in the figures. 

The fitted value of the restitution coefficient is e = 
0.87, in reasonable agreement with the approximate mea- 
surement described in Sec. Ill Bl The fitted values of the 
dimensionless heat flux at the sample bottom q*(0) are 
listed in Table [I] The total power input is computed 
as P — (ttD 2 /4)<Zz, where D — 0.87 cm is the effective 
sample diameter (Appendix . By equating the rate 
of momentum transfer to the total sample weight Nmg 
the input power is rigorously related to the average over 
grain impacts of the container-bottom velocity (wj) by 
P = 2Nmg(vi) 49]. Table[I]lists (vi) derived in this way 
from the fitted heat current, along with the ratio of (vi) 
to the peak container velocity vq = luzq. 

The fitted power input ql(0) is an increasing function 
of r, as expected on physical grounds. The power in- 
put also increases with bed depth Ni, which is expected 
from the momentum-balance argument, but the increase 
is very slow. One caveat is that these q* values are 
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FIG. 12: Density and temperature profiles for T = 18, with 
other conditions the same as in Fig. As in Fig. QE| the 
density at the bottom of the sample is much less than the 
peak density, and the measured granular temperature does 
not agree well with the theoretical granular temperature near 
the sample bottom. Note that this situation is achieved here 
by increasing V relative to the reference state, while in Fig. 1101 
it was achieved by decreasing Ne. 



TABLE I: Fitted dimensionless heat flux into the bottom 
of the sample g* (0) as a function of dimensionless vibration 
acceleration V and bed depth Ne. Also listed are the average 
container impact velocity (vi) derived from q* z (0) and the ratio 
of {vi) to the vibration velocity amplitude Vq. 
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(Vi)/v 


15 
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15 
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3.5 
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0.110 


18 


3.0 


4.8 


9.4 cm/s 


0.168 



derived from the hydrodynamic fits, which tend to fail 
where the density at the sample bottom is much less than 
the peak density fFigs. HUll2l) . In these cases particularly, 
the fitted q* represents an effective power input reflect- 
ing the bottom boundary condition, which is greater than 
the actual power input because the hydrodynamic fit has 
larger granular temperature than the experimental data 
and hence greater collision losses in the lower part of the 
sample. 
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One robust conclusion from Table [I] is that the av- 
erage container-grain impact velocity is much less than 
the peak container velocity, (Vi) <C vq (substituting the 
actual input power for the effective power only strength- 
ens this result). This conclusion corresponds to the ex- 
perimental observation that the sample "floats" above 
the vibrating container bottom (Fig. |2J, so that most 
container-grain impacts occur near the top of the vibra- 
tion stoke where the container velocity is much smaller 
than its peak value. 

In both the data and the hydrodynamic fits in Figs. El- 
El there is a "temperature inversion" — the temperature 
increases with height near the free upper surface of the 
sample. In the hydrodynamic theory the heat flux is 
strictly upwards so the temperature inversion is only 
made possible by the term in /i in Eq.0 The temperature 
inversion may be visible in some of the earliest experi- 
ments on quasi-two dimensional systems l20U2ll and has 
been extensively discussed theoretically £a> S Ha] The 
results presented here are perhaps the first demonstration 
that the temperature inversion in a physical system can 
be quantitatively explained by granular hydrodynamics. 



V. CONCLUSIONS 

One conclusion of this work is that it is possible, un- 
der certain conditions, to quantitatively model a physical 
vibrofluidized granular system using granular hydrody- 
namics. This is true despite the small sizes of the systems 
measured in grain diameters: the systems were small hor- 
izontally due to apparatus constraints, and in any case 
were necessarily small vertically in order to meet the flu- 
idization condition X < 1. The overall success of the 
hydrodynamic description, demonstrated by the fits in 
Figs. 191121 can probably be attributed to the fact that hy- 
drodynamics only breaks down when the mean free path 
becomes much larger than other relevant length scales. 

A limitation of present hydrodynamic theory is evi- 
dent in these figures — the inability to accurately describe 
the complex process of energy transfer from the vibrat- 
ing container bottom to the sample. The actual velocity 
distribution is highly skewed and anisotropic near the 
sample bottom (Fig.^J, and a discrepancy develops be- 
tween the theoretical and measured granular tempera- 
tures when the vibrational power input is large enough 
to deplete the particle density near the sample bottom 
(Figs HOI 113 ■ 

To calculate the power input, typically it has been as- 
sumed that the mean free path is much larger than the 
vibration stroke, or that the vibration waveform is saw- 
tooth or triangular |49l IHtI ] . Such waveforms are unreal- 
istic since they require infinite acceleration at the peak 
of the vibration stroke, which is precisely the point at 
which energy transfer to the grains occurs (Fig.|3J). Sim- 
ilarly the mean free path is much less than the vibration 
stroke in the current experiments (Appendix [5J|, which 
will be true in general for dense uniformly fluidized sys- 



tems unless very high vibration frequencies are used. One 
possibility for deriving the power input might be to use 
time- dependent hydrodynamics to model the formation 
and dissipation of shocks, as in Ref. l43l 

Thus it appears that current granular hydrodynamic 
theory can accurately describe the uniformly vibroflu- 
idized state well away from the vibrating energy source, 
but that derivation of the time-averaged boundary con- 
dition for realistic vibration waveforms remains an un- 
solved problem. If that problem can be solved, a 
complete first-principles description of the vibrofluidized 
state applicable to real granular systems will be in hand. 
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APPENDIX A: HORIZONTAL DENSITY 
PROFILES AND CALIBRATION OF BED 
HEIGHT 

The hydrodynamic model to which we have fit our data 
assumes that all quantities depend only upon the height 
coordinate z. This is not strictly true due to boundary 
effects at the vertical wall, as demonstrated by previous 
experiments |37J . The boundary condition also affects the 
effective cross sectional area of the sample tube, which 
is needed to compute the number of grain layers Ng for 
a given number of grains N. To address these issues 
we have taken some two-dimensional MRI data (Fig. I13|) 
to examine the horizontal density profile of the vibroflu- 
idized system, 

p(x,z) = J n(x,y,z)dy. (Al) 

If the NMR signal from each grain comes precisely from 
the grain center and the grain centers are uniformly dis- 
tributed horizontally over a cylinder of diameter D, the 
profile has the form 

p{XjZ)= im[l-(2 X /D)^ i {x <D/2, 
I otherwise. 

We also considered the possibility that the NMR signal 
from a grain comes from a sphere of diameter S < d, 
in which case the profile of Eq. IA2I must be convoluted 
with the profile for a single grain. We find that the fitted 
sample diameter is not significantly changed by varying 5 
over the entire possible range, Fig. 1131 In this figure some 
systematic deviation between the data and fit is visible, 
hinting that the density may be slightly depressed in the 
sample center 37] . However, these data are too noisy 
to be directly inverted to yield the radial density profile. 
The fitted diameter is consistent with a constant value 
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FIG. 13: Horizontal density profiles p(x,z) for z = 1.75 cm 
(top) and z = 1.55 cm (bottom), for the Nt ~ 3.0, T = 15 
reference state. The solid curves show fits to the profile for 
a uniform distribution of point-source grains in a cylindri- 
cal volume. The dashed curves show fits assuming the grain 
centers are uniformly distributed, but that the NMR signal 
source is a sphere with the full grain diameter. The differ- 
ence in the fitted cylinder diameter between the two fits is 
not significant. 



D = 8.7 mm (Fig. ll4Jl . This is less than the actual inside 
diameter (9 mm) but by less than the grain diameter 
(1.84 mm), suggesting a slight density enhancement at 
the walls due to the dynamic boundary condition. 



For the fits to the hydrodynamic theory, the height of 
the granular bed is specified by the dimcnsionlcss number 
density n* A — uaJ 2 — Nd 2 /(irD 2 A second measure 
of the bed depth is the number of layers Ne, defined as 
the height of the bed at rest divided by the grain diam- 
eter. The value of Nj> depends on the assumed volume 
filling fraction of the bed at rest (f>o via Ng = (ir /6(po)n\. 
For the Nt values quoted in this paper, shown in Ta- 
ble [HJ we have used 6p = 0.60 corresponding to a "loose 
random packing" [20|. It should be emphasized that the 
hydrodynamic fits are independent of the assumed value 
for cf>o, as they are based on n* A and not Ng. 




FIG. 14: Fitted effective tube diameter D as a function of 
height z for the Ni = 3.0, V = 15 reference state. The error 
bars become large at the top and bottom of the sample where 
the density is low. The dashed line shows the average value 
D = 8.7 mm used to analyze the data. 



TABLE II: Three measures of the sizes of the granular sam- 
ples used in this study. The number of grains TV was counted 
precisely. The dimensionless area density n* A is computed 
from N, the measured average grain diameter, and the fitted 
effective sample tube diameter from MRI experiments on the 
reference vibrofluidized state. The number of layers Nt de- 
pends additionally upon the assumed volume filling fraction 
of the bed at rest. The quantity X which characterizes the 
possibility of fully fluidizing the samples is also listed, using 
the value e = 0.87 resulting from the hydrodynamic fits. 
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Nt 


X = Ni(l-e) 
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3.42 


3.0 


0.39 
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4.56 


4.0 


0.52 



APPENDIX B: THE MEAN FREE PATH 

For hydrodynamics to be valid, the mean free path A 
must not be too large relative to distances over which 
the hydrodynamic fields vary. Therefore it is important 
to estimate A for the experimental conditions, and this 
also provides a consistency check for the measured bal- 
listic/diffusive crossover time. Extending the estimate of 
Ref. to three dimensions gives 



A* 



X/d 



A(n* 



n*{B - n*) 



(Bl) 



where A and B are chosen to give the exact dilute limit 
A* = l/(7T\/2w*) for n* — > and the heuristic dense limit 
A* = (1 - n*/n^)/3 for n* -> n%. Here n* Q = (6/tt)0 o is 
the density at which the grains come into contact, which 
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FIG. 15: Mean free path A as a function of number density 
n computed from Eq. IB1I Both quantities are nondimension- 
alized by the grain diameter d. 



we arbitrarily assume occurs at the volume filling fraction 



for loose random packing (f> Q = 0.60. 

The mean free path A*(n*) computed from Eg. IB II is 
shown in Fie. 1151 At the peak density n* = 0.85 that oc- 
curs in our experiments (Fig. Illjl . Ea. IB II gives A* 0.1. 
Thus it is hardly surprising that the hydrodynamic the- 
ory can explain the spatial dependence of the density and 
granular temperature in the high-density regions of the 
sample. Conversely, the data fit the hydrodynamic the- 
ory well for densities extending below n* — 0.1, which 
corresponds to A* > 2. Thus, the theory continues to 
work for mean free paths comparable to or even some- 
what greater than the scale over which the temperature 
and density change. 

The mean free path calculation can also be used to 
check the consistency of observed the ballistic crossover 
time, Fig. [5j For example, the center panel in this 
figure corresponds to an average dimensionless height 
z* = 5.2, and at this height the measured dimension- 
less density and temperature are n* = 0.41, T* = 0.5 
(Fig. EJ). Using A = d\*{n*) and v rms = (ST*gd)^ 2 
we calculate |6J| the ballistic to diffusive crossover time 
t CT = (32/37r) 1 / 2 A/?; rms = 6.3 ms, in excellent agreement 
with the observed value (Fig. |SJ . 
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